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Abstract 

We extend the one-body phase function upper bound on the superfluid fraction f s in a periodic 
solid (a spatially ordered supersolid) to include two-body phase correlations. The one-body current 
density is no longer proportional to the gradient of the one-body phase times the one-body density, 
but rather it becomes j(r{) = Pi(n)^\7 1 (p 1 (fi) + jf f cF 2/ o 2 (ri, f 2 )^Vi</> 2 (fi, f 2 ). This expression 
therefore depends also on two-body correlation functions. The equations that simultaneously 
determine the one-body and two-body phase functions require a knowledge of one-, two-, and 
three-body correlation functions. The approach can also be extended to disordered solids. Fluids, 
with two-body densities and two-body phase functions that are translationally invariant, cannot 
take advantage of this additional degree of freedom to lower their energy. 

PACS numbers: 67.80.-s, 67.90,+z, 67.57.De 
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I. INTRODUCTION 



id samples, but only in one case is it reported that the effect 
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In 1970, Leggett suggested that solid 4 He might exhibit non-classical rotational inertia 
(NCRI).[_lJ Recently, Kim and Chan have observed this effect for solid 4 He in Vycor, porous 
gold, and in bulk.[2j, y, |j] Given that, despite the extremely different porosities of these 
samples, the extrapolated values of the temperature T = K NCRI fraction (or NCRIF) 
are all in the 1-2% range, it seems unlikely that the effect occurs only at surfaces. Other 
experiments reproduce these results qualitatively, but with slightly smaller values, according 
to how the sample is prepared. In all previous experiments annealing has been used 

to improve the quality of the so 
can disappear under annealing. 

We have recently computed, using one-body densities obtained with an exact T=0 K 
Quantum Monte Carlo method, the superffuid fraction upper bound due to Leggett, based 
on a one-body phase function and the one-body density [81]. Our values were of the order of 
20%, far above what has been observed. If supersolidity indeed does occur, there is a need to 
understand why it does so at much lower values. For that reason, and because superfluidity 
in presence of a traslational broken symmetry may occur under other circumstances (such 
as condensates in optical lattices), we have re-examined the question of superfluidity in 
a supersolid, permitting the system to develop a phase function that includes not merely 
one-body effects, but also two-body (and, in principle, higher) terms. 

Leggett 's original study considered superflow in a solid confined to an annulus. He first 
ruled out systems whose ground state wavefunction ty is "disconnected" when a single 
particle coordinate is moved around the annulus (a generalized insulating state), and systems 
for which the wavefunction re-organizes to a lower energy state (a generalized form of the 
Landau criterion). To obtain an upper limit on the NCRIF, he then explicitly considered 
that, under rotation, the wavefunction became \I/ = ^oexp(i(j)), with a phase function cf) that 
was a sum of terms depending only upon a single coordinate (what we call a one-body phase 
function). pj He noted that this would yield only an upper limit to the T = K superffuid 
fraction f s . Slightly prior to Leggett, Andreev and Lifshitz posited superfluidity due to 
quantum diffusion of unspecified defect excitations. 9} A macroscopic theory of superflow 
was developed, but no microscopic definition of the superffuid velocity was given. At about 
the same time Chester [1Q] noted the existence of variational ground state wave functions 
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which display translational broken symmetry, corresponding to a crystal with vacancies; for 
such wave functions Bose Einstein condensation (BEC) had been demonstrated 0| . Kef.! 
suggested that capillary flow with a fixed lattice might occur under gravity, and one might 
expect other possibilities for unusual flow associated with solid 4 He. Despite numerous 
efforts, no such flow has been observed. 
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I3I [y], Q Additional, primarily topological, 



considerations have been given by Leggett. 



16] who in more recent work introduces the 



symbol f s for the superfluid fraction (SFF).[17| If superflow is the correct interpretation of 
the NCRI experiments, then f s is equivalent to the NCRIF. 

At the microscopic level, a number of Monte Carlo calculations for solid 4 He have been 
performed, with no common conclusion. Ceperley and Bernu have performed calculations 
on solid 4 He at low finite temperature with path integral Monte Carlo computations and 
to date they find evidence that is contrary to superfluidity. [18] They employ the winding 
number method, which yields reasonable values for the superfluid properties of liquid 4 He. 
For smaller samples they find a finite f 8 , but for larger samples they find f s = 0, although 
they observe that this negative result from the winding number calculation for the larger 
samples is inconclusive, due to possible phase-space sampling problems. They also calculate 
the exchange integrals, which their method successfully obtains for solid 3 He, and they 
interpret their results as arguing against superfluidity. Prokofiev and Svistunov also argue 



against bulk superfluidity. 



191 ] On the other hand, Galli, Rossi, and Reatto studied off- 



diagonal long-range order (ODLRO) at zero temperature for relatively large systems on the 
basis of an accurate variational theory, finding small but definite asymptotic values, which 
imply a finite BEC in the zero-momentum state, and "connectivity" of the wavefunction. {20I 
As a consequence, a superfluid velocity can be defined, and one can expect superflow if the 
excitations satisfy the Landau criterion. More recently, it has been found that superfluidity 
can occur in Monte Carlo calculations for a disordered "sample" . [2l| Such a sample was not 
equilibrated, and was deliberately disordered, and thus is far from equilibrium. One cannot 
rule out the possibility that this occurs experimentally, although it is not clear that the 
experimental samples are as disordered as the simulated samples. An exact computation 22j 
at zero temperature for a crystal with a finite concentration of vacancies gave ODLRO 
and therefore BEC but again it is not known if the experimental samples contain a finite 
concentration of vacancies at the low temperature of the experiment. 

Having ODLRO and satisfying the Landau criterion should ensure superfluidity, but it 
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still leaves open the problem of computing the value of f s ; this amounts to finding the 
best phase function, in the sense of minimizing the energy associated with superflow, where 
the phase gradient gives the superfluid velocity. The best such one-body phase function 
(but perhaps not the best phase function) satisfies the continuity equation with the current 
given by the product of the one-body density and the one-body superfluid velocity. [23| This 
provides an equation that enables the superfuid velocity profile to be determined from the 
average superfluid velocity and the local one-body density profile. Using a model for hep 
4 He as an fee lattice of the same density, ref. |23| took the one-body density to be a sum of 
gaussians. For relatively localized gaussians, the original calculations were not converged, 
due to computer limitations at that time, and thus provided only an upper limit for the 
superfluid fraction from the optimized one-body function. Using a 1976 estimate of the 
best gaussian, and more extensive computation, ref. |24j gave an upper bound for f s , or 



NCRIF, on the order of 2%, apparently in good agreement with experiment. However, 
more accurate Gaussian approximations to the one-body density, determined from "exact" 
Monte Carlo one-body densities, and the one-body densities themselves, have recently been 
shown to give f s on the order of 20% for solid 4 He at density of 0.029 A~ 3 (near the melting 
pressure) . |§j This is clearly too large a value, and requires re-examination of the theory. In 
what follows we therefore consider the possibility that a superfluid can respond to rotation 
by developing a phase that also includes two-body (and, in principle, higher) terms. In this 
case the superflow for each particle reflects not merely the local number density p(r), but 
also correlations with the other particles. 

Thus, as a particle flows around the crystal, it sees not only the average mass density, 
but also details of how other particles correlate to it, and changes its phase in response 
to them. When the two-body phase term is included, the one-body current density, which 
specifies the superflow pattern, is more complex than the usual j(f) = p{f)v s {f) form; it 
now depends explicitly on the two-body density (see Eq. [T2l) . and implicitly on both the 
two- and three-body densities. Although the former is occasionally computed (for example, 
in correlated basis-function studies), the three-body density is not normally studied. In the 
absence of such information we merely present the formalism, and do not perform actual 
calculations. Neverthess, this extra degree of freedom should enable the solid to lower its 
flow energy, and therefore to decrease the f s . This extension of the phase function to include 
two-body and higher terms should also apply to other candidates for supersolid behavior, 
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such as bosons confined in optical lattices. 

For bulk superfluid, the translational invariance of certain integrals involving the two- 
body correlation function and the two-body phase function prevents the system from utilizing 
the two-body phase term. This explains why, for a bulk superfluid at T = 0, theory based 
on a one-body phase alone gives the maximum possible value for f s , which is unity. 

Note that, if superfluidity is a bulk phenomenon in disordered Bose systems, then with 
certain developments the method of the present work should give the two-body upper bound 
for f s . The advantage of having a means to compute also upper bounds on f s , as opposed 
to only Monte Carlo calculations of f s , is that it ensures consistency and reveals the extent 
to which correlations are involved in the superfluidity. 



II. ENERGY OF SUPERFLOW FOR A NON-UNIFORM SYSTEM WITH N- 
BODY CORRELATIONS 

In what follows Xj will be used to represent the vector position x« of the i-th particle. 
Consider N identical Bose particles of mass m, with ground-state wavefunction ^o( x i • • • x n) 
and energy E Q . is an eigenstate of the Hamiltonian 

N 

2in 



where V is a sum of two-body interactions. Because \l/o is the ground state, it can be 
considered to be a real function, up to a constant phase. Our goal is to find the energy E 
of the system when there is an imposed average superflow velocity Vq. The energy E — E 
provides a measure of the superfluid fraction. 

Specifically, we consider superflow for a periodic geometry with average velocity Vq, with 
macroscopic repeat distance L and cross-section A. This imposes a phase change of mvoL/h 
on taking a single particle through L. Following Leggett, we assume that the effect on the 
wavefunction is simply to change its phase, via 

$ _> = ^ o exp(i0), (2) 

lion. Leggett 



where <p will depend on all N coordinates, perhaps in a rather complex fas 
has also considered additional possibilities that we do not discuss here. [16 



17) Note that 
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for such a response neither the density nor the potential energy changes; only the kinetic 
energy changes. 

With normalization f \ty Q \ 2 dxi . . . dx^ = 1, it is straightforward to show that 

K 2 N r 

AE = E-E = —J2j Itfol 2 ^) 2 ^! . • • dx N . (3) 

i 

Minimization of AE with respect to a that depends on all the particle coordinates (iV-body 
correlations) leads to Vj • (|\I/o| v »0) = for each i. 



III. SUPERFLOW FOR A NON-UNIFORM SYSTEM WITH TWO-BODY COR- 
RELATIONS 



Now consider a that depends on correlations between the coordinates. Let us expand 
it in the form 

N I N 

<P = Y1 fofc) + 2 (N - 1) ^ + "' ^ 

i i^j 

The one-body term includes both the imposed flow and the one-body backflow. The two- 
body contributions to the phase represent correlated backflow, and dots represent more 
complex backflow effects. By symmetry, <p2{xi,Xj) = 4>2{xj,Xi). There is no constant term 
in either coordinate, since that would be part of <p\. Hence, in the fourier expansion of 02, 
there are no zero-wavevector terms. Many-body phases have been considered previously, 
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27 



28 



but not in the context of steady bulk superflow. 

In terms of 0, which must be symmetric in all coordinates, we have 

h 1 N 

v sl (x u ...) = -V x = v s (l) + - — - V v s (l, i) + . . . , (5) 
m JM — 1 

¥i 

where v s {\) = ^Vi0i(xi) and v s (l, i) = ^Vi02(^i, %i) are expected to be of the same order 
of magnitude. If the N — 1 two-body velocities are to yield a net velocity that is comparable 
to those from the one-body velocity, then the two-body velocity should be on the order of 
the one-body velocity. 
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We then have 

N , N 



v 2 



1 N N 

The first term is what one would obtain from a one-body phase. The second is a (new) 
cross-term, involving both the one-body and two-body phases, that makes the one- and 
two-body phases dependent on one another. (There are N — l equal terms of this sort.) The 
third is the square of the two-body term involving 1 and i. (There are N — l equal terms 
of this sort; in the limit N — > oo, this term is negligible so that it is dropped from now on.) 
The fourth term gives the two-body phase cross-correlations. (There are 2(iV — 1)(N — 2) 
equal terms of this sort, so it dominates the contribution of the third term.) 

It is now convenient to define a number of quantitites. The one-body density p\ is given 

by 

pi(l) = pi(xi) =N J \^ Q \ 2 dx 2 . . . dx N , (7) 

where f dx\pi(x\) = N. The two-body density p 2 is given by 

p 2 (l, 2) = p 2 (x 1 ,x 2 ) = N(N -I) J |^o| 2 ^ 3 • • • dx N , (8) 

where f p 2 (xi,x 2 )dx 2 = (N — l)p 1 (x 1 ). The three-body density p 3 is given by 

p 3 (l,2,3) =p 3 ( Xl ,x 2 ,x 3 ) = N(N-l)(N-2) J \^ \ 2 dx 4 . . . dx N , (9) 

where J p 3 (xi,x 2 ,x 3 )dx 3 = (N — 2)p 2 (xi,x 2 ). With these conventional definitions, one has 
PnO-i ■ ■ ■ ,N) = AM|\I/o| 2 , which leads to additional factors of order N that later appear in 
the theory of the flow energy. For n <^ N, note that p n ~ (N/V) n . 

In terms of these densities, and in the limit iV — > oo, if the one-body and two-body 
velocities are of the same order of magnitude, we obtain for the flow energy 

AE -»■ — / p^x^Vs^fdxi + — / p 2 (x!, x 2 )v s (l) ■ v S)2 dx x dx 2 

Tft f 

+ 2iV2 / '° 3 ( a:i ' X2 ' X3 )^. 2 ' v s , 3 dxidx 2 dx 3 . (10) 
The second and the third terms of this equation represent correlated backfow contributions 
the flow energy. In the many-body expansion of 4> we have truncated correlated backflow to 
the pair level; these results would be more complex if triplet correlations were included in 
0. 



IV. MINIMIZATION OF THE ENERGY OF SUPERFLOW 



The condition that E be a minimum relative to variations in 4>i(xi) leads to 

0= h~5<jh = ^ 1 '^ Xl ^ s ^ + n J dx 2p2(x 1 ,x 2 )v s (l,2)} 

= Vi-/(xO, (11) 

— * 

where the one body current density j(xi) now is the sum of the usual one-body term and 
an integral over the two-body density: 

j(xi) = pi(xi)v 8 (l) + J dx 2 p 2 (x 1 ,x 2 )v s (l,2). (12) 

The condition that E be a minimum relative to variations in (f) 2 (xi,x 2 ) leads to (in the 
limit N — > 00) 

= YJ^~ = Vi • [P2(xi,x 2 )v s (l) + J dx 3 p 3 (x 1 ,x 2 ,x 3 )v s (l,3)]. (13) 

With V the system volume, and with periodic boundary conditions we expand the known 
densities via 

Pi( x i) = ^2p\cxp[iG x ■ £1], (14) 

A 

N ^ 

p 2 (x 1 ,x 2 ) = — P\,n ex P[^(^A ' X\ + G M • x 2 )}, (15) 
N 2 

^3(^1,2:2,^3) = ^ ^ P\,n,v exp[i(Gx ■x 1 +G l ^-x 2 + G u - x 3 )] . (16) 

Xfiu 

If the system is in a crystalline phase and the center of mass is fixed, then the densities 
must reflect the periodicity of the crystal lattice. In this case the vectors G\, G\ + G M , and 
G\ + G^ + G v must be vectors of the reciprocal lattice. These definitions give all the fourier 
coefficients (e.g. px, Px,^) the same dimensionality. Note that p = N/V is the average 
particle density. 

We now expand the unknown phases via 

fafa) = mV ° Xl + ^2<p\exp[iG x • xi], (17) 

n A^O 

(p 2 (x 1 ,x 2 ) = ^2 0A, M exp[i(G A • xi + G M • f 2 )]- (18) 
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— » 

Alternatively, in (I17p we can sum over all A by taking Go to be very small and 0o to be very 
large, subject to the condition that v s o = (ih/m)Go(fio be finite. The first term in (TTTT) is the 
imposed flow, and the second term is the (to-be-determined) one-body backflow. 

In practice, we consider v to be known, and we must find the other components in terms 
of vq or, equivalently, 0o- The unknown velocities become, via v s = (h/m)V(f), 

v s (xi) = vq + — "S^ G X (j) X exp[iG x ■ x x \ 
m 

A^O 

= v + ^2v x exp{iGx- xi], (19) 

A^O 

v s (x 1 ,x 2 ) = — G X (px^exp[i(G x ■ x\ + ■ x 2 )\ 
m ^— ' 

= E v x ^exp{i(Gx • x\ + G M ■ x 2 )), (20) 

A^0,^0 

where v x = {ih/m)G x <f) X , and v x ^ = [ih/m)G X (j) X ^. 
Eq. ( TTTj) implies that, for all A ^ 0, that 

(21) 

For a discretization where A takes on P values, this gives P — 1 conditions on P — 1 unknown 
values of (fi x (where 6q is considered to be known). For 0^,-7 = 0, these are the equations 



employed in refs. 



23 



24 



25 



Eq. (fl3|) implies that, for all A 7^ and all v ^ 0, that 

= E(^a • d»)\px-w<j>v + Epa-m.^-t]- (22) 

Note that all of the 0's are driven by the superflow, and thus must be proportional to Vq 
dotted with a vector associated with the system. For a discretization where A takes on P 
values and v takes on P values, this gives (P — l) 2 conditions on [P — l) 2 unknown values 
of M _ 7 . 
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V. SUPERFLUID DENSITY 



We now derive an expression for the superfTuid density. (Recall that this is really an 
upper limit to the superfluid density) Let us rewrite (JTOj) . for N — ► oo, as 

in f I / 

AE = — / dxiv s (l) ■ \pi (xi)v a (l) + — / dx 2 p 2 (xi,x 2 )v s (l,2)} 

Tfl f 

+ 2N / dxidx2 ^ s ( 1 ^ 2 ) ' v s (l)p2(xi,x 2 ) 

Tfl f 

+ 2iV2 / ^1^2^3^(1,2) • v s (l,3)p 3 (x 1 ,x 2 ,x 3 ). (23) 

Because v s (l) and v s (l,2) involve Vi, in both of these terms we can integrate by parts, 
where we take the integration direction along the flow. Each volume term then leads to 
a surface term and to a volume term that is the product of a phase with a divergence. 
However, by (jlip and ffTSl . these divergences are zero. Moreover, the quantity v s (l,2) is 
strictly periodic, so its phase is periodic, and therefore its two surface terms cancel. All that 
remains is the surface term associated with the change in phase of v s (l) as a particle goes 
around the system. Defining Si to be this surface, of area A, and defining vq to be the unit 
vector normal to Si, i.e. the direction of the imposed velocity ■uo, we then have 

AE = 7; J dAx(A(j)o)vo • [pi(ziH(l) + J dx 2 p2(x 1 ,x 2 )v s (l,2)] 

= ^A(A<j) )j = ^-ALv j = ^-V Ps vl (24) 

where V = AL and jo = Ps^o is given in fourier component form as 

= Povo + ^^G M [p_ M /t + P-,xjL<f>^-i\- (25) 

ft^O 1^0 

There is no contribution to jo from the part of ^2(^1,^2) that is independent of x± because 
the integral over x\ involves a strictly periodic function of x\. Hence, once we know the flow 
pattern in terms of v , we can compute j and determine p s = jo/v . The first term in (f2"5|) 
is the imposed flow, the second term is the one-body backflow, and the third term is the 
two-body backflow. 

The theory is suflicently general that it can be applied both to liquid and solid boson 
systems. For a homogeneous liquid the previous equations simplify drastically. In the liquid 
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the two-body density p 2 (xi, ^2) is a function only of \x\ — x 2 \- For the liquid we expect 
the two-body phase 4>2{x\,x 2 ) also to be a function only of \x\ — x 2 \. As a consequence, 
^2(^1,^2) is even in x\ — x 2 , but v s (l,2), involving Vi|a?i — x 2 \, is odd in x\ — x 2 . Hence 
the second term in (|24[) goes to zero. Therefore the two body terms do not enter explicitly 
in determining p s . Moreover, in fill I) , the two-body terms similarly are zero, so that they 
do not even affect the one-body terms. Hence the two-body terms, i.e. backflow, are 
irrelevant both to determining the one-body flow pattern and the total energy of a uniformly 
flowing translationally invariant system. Therefore it is valid to neglect two-body terms in 
determining the superfluid density in a fluid, as is done conventionally. We believe that 
study of the higher-order phase correlations will reveal that none of them are relevant to a 
superfluid. The one-body phase alone then gives a T = K superfluid fraction of unity, in 
agreement with a conclusion in ref . [17 . 



The actual computation of the upper bound on f s which includes the two-particle phase 
requires the solution of eq.(29), where the phases M and </> M ,- 7 are solutions of eqs.(25) and 
(26). As inputs one needs to know p\, p\ tfJl , and px,n,u, the fourier components of the one- 
body, two-body, and three-body correlation functions, respectively, as defined in eqs. (18-20). 
There is no difficulty in obtaining highly accurate values for p\ from an exact computation 8j 
for the ground state of solid 4 He. Also, the computation of p\ ttl = < Pq x Pq > is currently 
performed in Monte Carlo simulations and there is no difficulty in computing p\ tfl for the 
first few stars of the reciprocal wave vectors. The most complex terms are the three-body 
terms p\^ v . In principle, the computation of p\^ v is also straightforward with Monte Carlo 
methods, but in this case problems might arise in terms of the accuracy of such three-body 
correlation functions and the number of terms necessary in order to get convergence for the 
bound on f s . Only an explicit computation will be able to answer this. Often one avoids the 
explicit computation of three-body correlations by a suitable approximation in terms of cor- 
relation functions of lower order - for instance, the so-called convolution approximation. 29] . 
In conclusion, with present computational resources it should be possible to perform the 
computations needed to get an improved bound on / s . 
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VI. CONCLUSIONS 



Ref. [8| calculated an upper bound for the superfluid fraction using an optimized one-body 
phase function and state-of-the-art Monte Carlo one-body densities, finding an upper bound 
for the superfluid fraction of 20% for solid 4 He near the melting pressure, significantly higher 
than has been observed. Given an additional degree of freedom by permitting a many-body 
phase function, it should be possible to lower the flow energy for the case of the solid, and 
therefore to lower the corresponding superfluid fraction, perhaps by a considerable amount. 
The present work develops the theory for the superfluid fraction in a periodic solid including 
a two-body phase function. The results should also be applicable to condensates in optical 
lattices, and the approach should also apply to disordered solids. We have shown which 
quantities are needed for such a computation, which should give information on backflow 
effects in the solid phase. Computation of the necessary inputs - the fourier components 
of two- and three-body correlation functions - appears feasible with present computational 
resources, but this is left as a future development. 
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